#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
SCALE GROUP: COMPLETE NUMERICAL VERIFICATION SUITE
Formal verification of all theorems and properties

Author: Raheb Ali Mohammed Saleh Aoudh
Date: February 2026
Version: 3.0 (Final Version with Corrected Roots)

This code provides comprehensive numerical verification for the paper
"The Scale Group: A Novel Abelian Group Structure on the Positive Reals".
All theorems and identities are tested with high precision.
The root function has been corrected using the proper mathematical derivation.
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import special
import warnings
warnings.filterwarnings('ignore')

# =============================================================================
# CONFIGURATION
# =============================================================================

class Config:
    """Configuration parameters for numerical verification"""
    KAPPA = 10.0  # Default κ value (any positive number works)
    N_TESTS = 1000  # Number of random tests
    N_ZETA_TERMS = 10000  # Terms for zeta series
    PRECISION = 1e-12  # Precision threshold
    
config = Config()

# =============================================================================
# ERROR ANALYSIS CLASS
# =============================================================================

class ErrorAnalysis:
    """Comprehensive error analysis utilities"""
    
    @staticmethod
    def absolute_error(left, right):
        """Absolute error: |left - right|"""
        return abs(left - right)
    
    @staticmethod
    def relative_error(left, right, floor=1.0):
        """Relative error: |left - right| / max(|left|, |right|, floor)"""
        denominator = max(abs(left), abs(right), floor)
        return abs(left - right) / denominator if denominator > 0 else float('nan')
    
    @staticmethod
    def summary(errors):
        """Statistical summary of errors"""
        if not errors or all(np.isnan(errors)):
            return {
                'mean': float('nan'),
                'max': float('nan'),
                'min': float('nan'),
                'std': float('nan'),
                'count': 0
            }
        # Filter out nan values
        errors_clean = [e for e in errors if not np.isnan(e)]
        if not errors_clean:
            return {
                'mean': float('nan'),
                'max': float('nan'),
                'min': float('nan'),
                'std': float('nan'),
                'count': 0
            }
        return {
            'mean': np.mean(errors_clean),
            'max': np.max(errors_clean),
            'min': np.min(errors_clean),
            'std': np.std(errors_clean),
            'count': len(errors_clean)
        }
    
    @staticmethod
    def print_summary(title, errors):
        """Print formatted error summary"""
        print(f"\n  {title}:")
        if errors['count'] == 0:
            print(f"    No valid tests")
            return
        print(f"    Mean error: {errors['mean']:.2e}")
        print(f"    Max error:  {errors['max']:.2e}")
        print(f"    Min error:  {errors['min']:.2e}")
        print(f"    Std dev:    {errors['std']:.2e}")
        print(f"    Count:      {errors['count']}")

errors = ErrorAnalysis()

# =============================================================================
# SCALE GROUP CLASS (COMPLETELY FIXED)
# =============================================================================

class ScaleGroup:
    """
    Implementation of the scale group (M_κ, ⊗_κ)
    All properties are proved for any κ > 0
    """
    
    def __init__(self, kappa=config.KAPPA):
        """Initialize scale group with parameter κ > 0"""
        if kappa <= 0:
            raise ValueError(f"κ must be positive, got {kappa}")
        self.kappa = kappa
        self.identity = np.exp(1/kappa)
        self.name = f"Scale Group (κ={kappa:.3f})"
    
    def otimes(self, a, b):
        """
        Scale multiplication: a ⊗ b = exp(κ ln a ln b)
        """
        if a <= 0 or b <= 0:
            return float('nan')
        
        log_a = np.log(a)
        log_b = np.log(b)
        exponent = self.kappa * log_a * log_b
        
        # Handle potential overflow/underflow
        if exponent > 700:
            return np.inf
        if exponent < -700:
            return 0.0
            
        return np.exp(exponent)
    
    def inverse(self, a):
        """
        Group inverse: a^{-1} = exp(1/(κ² ln a))
        """
        if a <= 0:
            return float('nan')
        if abs(a - 1.0) < 1e-10:
            return 1.0
        
        log_a = np.log(a)
        exponent = 1.0 / (self.kappa**2 * log_a)
        
        if exponent > 700:
            return np.inf
        if exponent < -700:
            return 0.0
            
        return np.exp(exponent)
    
    def T_transform(self, x):
        """
        T_κ(x) = ln(κ ln x) for x > 1
        """
        if x <= 1:
            return float('nan')
        
        inner = self.kappa * np.log(x)
        if inner <= 0:
            return float('nan')
            
        return np.log(inner)
    
    def T_inverse(self, y):
        """
        Inverse of T_κ: T^{-1}(y) = exp(e^y / κ)
        """
        exponent = np.exp(y) / self.kappa
        
        if exponent > 700:
            return np.inf
            
        return np.exp(exponent)
    
    def power(self, a, n):
        """
        a^{⊗ n} = exp(κ^{n-1} (ln a)^n)
        """
        if a <= 1:
            return float('nan')
        if n <= 0:
            return float('nan')
        
        log_a = np.log(a)
        exponent = self.kappa**(n-1) * (log_a)**n
        
        if exponent > 700:
            return np.inf
        if exponent < -700:
            return 0.0
            
        return np.exp(exponent)
    
    def power_repeated(self, a, n):
        """
        Compute a^{⊗ n} by repeated multiplication (for verification)
        """
        if a <= 1 or n <= 0:
            return float('nan')
        
        result = a
        for _ in range(1, n):
            result = self.otimes(result, a)
            if not np.isfinite(result):
                return np.inf
        return result
    
    def root(self, a, n):
        """
        a^{⊗ 1/n} = exp( (κ ln a)^{1/n} / κ )
        
        Derived from: T(a^{⊗ 1/n}) = (1/n) T(a)
        where T(x) = ln(κ ln x)
        
        Proof:
        T(a^{⊗ 1/n}) = (1/n) T(a)
        ln(κ ln(a^{⊗ 1/n})) = (1/n) ln(κ ln a)
        κ ln(a^{⊗ 1/n}) = (κ ln a)^{1/n}
        ln(a^{⊗ 1/n}) = (κ ln a)^{1/n} / κ
        a^{⊗ 1/n} = exp( (κ ln a)^{1/n} / κ )
        """
        if a <= 1:
            return float('nan')
        if n <= 0:
            return float('nan')
        
        log_a = np.log(a)
        
        # (κ ln a)^{1/n}
        inner = (self.kappa * log_a) ** (1.0/n)
        
        # (κ ln a)^{1/n} / κ
        exponent = inner / self.kappa
        
        if exponent > 700:
            return np.inf
        if exponent < -700:
            return 0.0
            
        return np.exp(exponent)
    
    def root_verification(self, a, n):
        """
        Verify that (a^{⊗ 1/n})^{⊗ n} = a
        """
        root_val = self.root(a, n)
        if not np.isfinite(root_val):
            return float('nan')
        power_check = self.power(root_val, n)
        return power_check
    
    def zeta_scale(self, s, terms=config.N_ZETA_TERMS):
        """
        ζ_κ(s) = Σ_{n=1}^∞ n^{-⊗ s} = Σ n^{-κ ln s}
        """
        if s <= 1:
            return float('nan')
        
        exponent = -self.kappa * np.log(s)
        result = 0.0
        
        for n in range(1, terms + 1):
            try:
                term = n ** exponent
                if not np.isfinite(term):
                    break
                result += term
                if term < 1e-15 and n > 100:
                    break
            except:
                continue
                
        return result
    
    def zeta_classical(self, s):
        """ζ(κ ln s) using scipy"""
        try:
            arg = self.kappa * np.log(s)
            if arg <= 1:
                return float('nan')
            return special.zeta(arg, 1)
        except:
            return float('nan')
    
    def scale_prime(self, q):
        """Scale prime corresponding to ordinary prime q: p = exp(e^q / κ)"""
        try:
            exponent = np.exp(q) / self.kappa
            if exponent > 700:
                return np.inf
            return np.exp(exponent)
        except:
            return np.inf


# =============================================================================
# 1. GROUP AXIOMS VERIFICATION
# =============================================================================

def verify_group_axioms(sg, n_tests=config.N_TESTS):
    """Verify all group axioms with random testing"""
    print("\n" + "=" * 80)
    print("1. GROUP AXIOMS VERIFICATION")
    print("=" * 80)
    
    # Generate random test points
    a = np.random.uniform(1.1, 50, n_tests)
    b = np.random.uniform(1.1, 50, n_tests)
    c = np.random.uniform(1.1, 50, n_tests)
    
    # 1.1 Closure
    print("\n1.1 Closure: a ⊗ b > 0")
    closure_ok = True
    for i in range(min(100, n_tests)):
        result = sg.otimes(a[i], b[i])
        if not (result > 0 and np.isfinite(result)):
            closure_ok = False
    print(f"  ✓ Closure: {closure_ok}")
    
    # 1.2 Associativity
    assoc_errors = []
    for i in range(n_tests):
        try:
            left = sg.otimes(sg.otimes(a[i], b[i]), c[i])
            right = sg.otimes(a[i], sg.otimes(b[i], c[i]))
            if np.isfinite(left) and np.isfinite(right):
                assoc_errors.append(errors.relative_error(left, right))
        except:
            continue
    assoc_summary = errors.summary(assoc_errors)
    errors.print_summary("Associativity", assoc_summary)
    
    # 1.3 Commutativity
    comm_errors = []
    for i in range(n_tests):
        try:
            left = sg.otimes(a[i], b[i])
            right = sg.otimes(b[i], a[i])
            if np.isfinite(left) and np.isfinite(right):
                comm_errors.append(errors.relative_error(left, right))
        except:
            continue
    comm_summary = errors.summary(comm_errors)
    errors.print_summary("Commutativity", comm_summary)
    
    # 1.4 Identity
    id_errors = []
    for i in range(n_tests):
        try:
            result = sg.otimes(a[i], sg.identity)
            if np.isfinite(result):
                id_errors.append(errors.relative_error(result, a[i]))
        except:
            continue
    id_summary = errors.summary(id_errors)
    errors.print_summary(f"Identity (e={sg.identity:.6f})", id_summary)
    
    # 1.5 Inverse
    inv_errors = []
    for i in range(n_tests):
        try:
            inv = sg.inverse(a[i])
            result = sg.otimes(a[i], inv)
            if np.isfinite(result):
                inv_errors.append(errors.relative_error(result, sg.identity))
        except:
            continue
    inv_summary = errors.summary(inv_errors)
    errors.print_summary("Inverse", inv_summary)
    
    return {
        'associativity': assoc_summary,
        'commutativity': comm_summary,
        'identity': id_summary,
        'inverse': inv_summary
    }


# =============================================================================
# 2. ISOMORPHISM VERIFICATION
# =============================================================================

def verify_isomorphism(sg, n_tests=config.N_TESTS):
    """Verify the isomorphism T: (M,⊗) → (ℝ,+)"""
    print("\n" + "=" * 80)
    print("2. ISOMORPHISM VERIFICATION")
    print("=" * 80)
    
    a = np.random.uniform(1.1, 50, n_tests)
    b = np.random.uniform(1.1, 50, n_tests)
    
    # 2.1 Homomorphism property
    print("\n2.1 Homomorphism property: T(a ⊗ b) = T(a) + T(b)")
    homo_errors = []
    for i in range(n_tests):
        try:
            ab = sg.otimes(a[i], b[i])
            if np.isfinite(ab) and ab > 1:
                left = sg.T_transform(ab)
                right = sg.T_transform(a[i]) + sg.T_transform(b[i])
                if np.isfinite(left) and np.isfinite(right):
                    homo_errors.append(errors.absolute_error(left, right))
        except:
            continue
    homo_summary = errors.summary(homo_errors)
    errors.print_summary("Homomorphism", homo_summary)
    
    # 2.2 Injectivity
    print("\n2.2 Injectivity: T(a) = T(b) ⇒ a = b")
    inj_errors = []
    for i in range(100):
        try:
            Ta = sg.T_transform(a[i])
            if np.isfinite(Ta):
                b_recovered = sg.T_inverse(Ta)
                if np.isfinite(b_recovered):
                    inj_errors.append(errors.relative_error(b_recovered, a[i]))
        except:
            continue
    inj_summary = errors.summary(inj_errors)
    errors.print_summary("Injectivity", inj_summary)
    
    # 2.3 Surjectivity
    print("\n2.3 Surjectivity: ∀ y ∃ x with T(x) = y")
    y = np.random.uniform(-5, 5, 100)
    surj_errors = []
    for i in range(100):
        try:
            x = sg.T_inverse(y[i])
            if np.isfinite(x) and x > 1:
                y_recovered = sg.T_transform(x)
                if np.isfinite(y_recovered):
                    surj_errors.append(errors.absolute_error(y_recovered, y[i]))
        except:
            continue
    surj_summary = errors.summary(surj_errors)
    errors.print_summary("Surjectivity", surj_summary)
    
    # 2.4 Inverse property
    print("\n2.4 Inverse property: T^{-1}(T(x)) = x")
    inv_prop_errors = []
    for i in range(100):
        try:
            Tx = sg.T_transform(a[i])
            if np.isfinite(Tx):
                x_recovered = sg.T_inverse(Tx)
                if np.isfinite(x_recovered):
                    inv_prop_errors.append(errors.relative_error(x_recovered, a[i]))
        except:
            continue
    inv_prop_summary = errors.summary(inv_prop_errors)
    errors.print_summary("Inverse property", inv_prop_summary)
    
    return {
        'homomorphism': homo_summary,
        'injectivity': inj_summary,
        'surjectivity': surj_summary,
        'inverse': inv_prop_summary
    }


# =============================================================================
# 3. POWERS AND ROOTS VERIFICATION (COMPLETELY FIXED)
# =============================================================================

def verify_powers_and_roots(sg, n_tests=50):
    """Verify power and root formulas"""
    print("\n" + "=" * 80)
    print("3. POWERS AND ROOTS VERIFICATION")
    print("=" * 80)
    
    a = np.random.uniform(1.1, 10, n_tests)
    exponents = [1, 2, 3, 4, 5]
    
    # 3.1 Powers
    print("\n3.1 Powers: a^{⊗ n}")
    power_results = {}
    for n in exponents:
        errors_list = []
        for ai in a:
            try:
                power_formula = sg.power(ai, n)
                power_repeated = sg.power_repeated(ai, n)
                if np.isfinite(power_formula) and np.isfinite(power_repeated):
                    errors_list.append(errors.relative_error(power_formula, power_repeated))
            except:
                continue
        power_results[f'n={n}'] = errors.summary(errors_list)
        errors.print_summary(f"Power n={n}", power_results[f'n={n}'])
    
    # 3.2 Roots (COMPLETELY FIXED)
    print("\n3.2 Roots: a^{⊗ 1/n}")
    print("    Using formula: a^{⊗ 1/n} = exp( (κ ln a)^{1/n} / κ )")
    root_results = {}
    for n in exponents:
        errors_list = []
        for ai in a:
            try:
                # Compute root using corrected formula
                root_val = sg.root(ai, n)
                if np.isfinite(root_val):
                    # Verify that (root)^{⊗ n} = a
                    power_check = sg.power(root_val, n)
                    if np.isfinite(power_check):
                        rel_error = errors.relative_error(power_check, ai)
                        errors_list.append(rel_error)
                        
                        # Print first few for manual verification
                        if len(errors_list) <= 3 and n == 2:
                            print(f"      a={ai:.2f}, root={root_val:.6f}, "
                                  f"(root)^n={power_check:.6f}, error={rel_error:.2e}")
            except:
                continue
        root_results[f'n={n}'] = errors.summary(errors_list)
        errors.print_summary(f"Root n={n}", root_results[f'n={n}'])
    
    return {'powers': power_results, 'roots': root_results}


# =============================================================================
# 4. ZETA FUNCTIONS VERIFICATION
# =============================================================================

def verify_zeta_functions(sg):
    """Verify scale zeta function properties"""
    print("\n" + "=" * 80)
    print("4. ZETA FUNCTIONS VERIFICATION")
    print("=" * 80)
    
    # 4.1 Fundamental identity
    print("\n4.1 Fundamental identity: ζ_κ(s) = ζ(κ ln s)")
    s_values = [1.2, 1.5, 2.0, 2.5, 3.0]
    
    print("\n  s     κ ln s    ζ_κ(s)      ζ(κ ln s)    error")
    print("  " + "-" * 60)
    
    zeta_errors = []
    for s in s_values:
        try:
            zeta_s = sg.zeta_scale(s)
            zeta_classical = sg.zeta_classical(s)
            
            if np.isfinite(zeta_s) and np.isfinite(zeta_classical):
                abs_error = errors.absolute_error(zeta_s, zeta_classical)
                rel_error = errors.relative_error(zeta_s, zeta_classical)
                zeta_errors.append(rel_error)
                print(f"  {s:.1f}   {sg.kappa * np.log(s):.2f}   {zeta_s:.8f}   {zeta_classical:.8f}   {abs_error:.2e}")
            else:
                print(f"  {s:.1f}   {sg.kappa * np.log(s):.2f}   {'NaN':<10}   {'NaN':<10}   NaN")
        except Exception as e:
            print(f"  {s:.1f}   {sg.kappa * np.log(s):.2f}   Error     Error     Error")
    
    # 4.2 Zeros
    print("\n4.2 Zeros of ζ_κ(s)")
    riemann_zeros_t = [14.1347, 21.0220, 25.0109, 30.4249, 32.9351]
    expected_radius = np.exp(1/(2*sg.kappa))
    print(f"\n  Expected radius: |s| = e^(1/(2κ)) = {expected_radius:.6f}")
    
    zero_errors = []
    for t in riemann_zeros_t:
        s = np.exp((0.5 + 1j * t) / sg.kappa)
        radius = abs(s)
        error = abs(radius - expected_radius)
        zero_errors.append(error)
        print(f"  t={t:.4f} → |s|={radius:.6f}, error={error:.2e}")
    
    zero_summary = errors.summary(zero_errors)
    print(f"\n  Mean error: {zero_summary['mean']:.2e}")
    print(f"  Max error:  {zero_summary['max']:.2e}")
    
    return {
        'zeta_identity': errors.summary(zeta_errors),
        'zeta_zeros': zero_summary
    }


# =============================================================================
# 5. SCALE PRIMES VERIFICATION
# =============================================================================

def verify_scale_primes(sg):
    """Verify scale prime properties"""
    print("\n" + "=" * 80)
    print("5. SCALE PRIMES VERIFICATION")
    print("=" * 80)
    
    # 5.1 Prime correspondence
    print("\n5.1 Prime correspondence: p = exp(e^q/κ), T(p) = q")
    primes = [2, 3, 5, 7, 11, 13, 17]
    
    print("\n  q    scale prime p           T(p)        error")
    print("  " + "-" * 60)
    
    prime_errors = []
    for q in primes:
        p = sg.scale_prime(q)
        if np.isfinite(p) and p != np.inf and p > 1:
            T_p = sg.T_transform(p)
            if np.isfinite(T_p):
                error = abs(T_p - q)
                prime_errors.append(error)
                print(f"  {q:2d}    {p:.6e}    {T_p:.6f}    {error:.2e}")
            else:
                print(f"  {q:2d}    {p:.6e}    {'NaN':<10}    NaN")
        else:
            print(f"  {q:2d}    {'∞':<15}    {'∞':<10}    ∞")
    
    # 5.2 Prime counting
    print("\n5.2 Scale prime counting: π_κ(x) = π(T(x))")
    
    def is_prime(n):
        if n < 2:
            return False
        if n == 2:
            return True
        if n % 2 == 0:
            return False
        for i in range(3, int(np.sqrt(n)) + 1, 2):
            if n % i == 0:
                return False
        return True
    
    def pi_ordinary(y):
        if y < 2:
            return 0
        return sum(1 for n in range(2, int(y) + 1) if is_prime(n))
    
    x_values = [2, 5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000]
    
    print("\n  x          T(x)     π_κ(x)   T/ln T     ratio")
    print("  " + "-" * 55)
    
    counting_results = []
    for x in x_values:
        if x > 1:
            Tx = sg.T_transform(x)
            if np.isfinite(Tx) and Tx > 1:
                pi_scale = pi_ordinary(Tx)
                asymptotic = Tx / np.log(Tx) if Tx > 1 else 0
                ratio = pi_scale / asymptotic if asymptotic > 0 else 0
                counting_results.append({
                    'x': x, 'Tx': Tx, 'pi_scale': pi_scale, 
                    'asymptotic': asymptotic, 'ratio': ratio
                })
                print(f"  {x:.2e}  {Tx:.2f}   {pi_scale:6d}   {asymptotic:.2f}   {ratio:.3f}")
    
    return {
        'prime_errors': errors.summary(prime_errors),
        'counting': counting_results
    }


# =============================================================================
# 6. ADDITIONAL TESTS
# =============================================================================

def verify_additional_tests(sg):
    """Run additional numerical stability and edge case tests"""
    print("\n" + "=" * 80)
    print("6. ADDITIONAL TESTS")
    print("=" * 80)
    
    # 6.1 Different κ values
    print("\n6.1 Testing with different κ values")
    kappa_values = [0.1, 1.0, 10.0, 100.0, 1000.0]
    a, b = 2.0, 3.0
    
    for kappa in kappa_values:
        try:
            sg_test = ScaleGroup(kappa)
            ab = sg_test.otimes(a, b)
            if np.isfinite(ab) and ab > 1:
                T_ab = sg_test.T_transform(ab)
                T_sum = sg_test.T_transform(a) + sg_test.T_transform(b)
                error = abs(T_ab - T_sum)
                print(f"  κ={kappa:6.1f}: T(a⊗b) - (T(a)+T(b)) = {error:.2e}")
            else:
                print(f"  κ={kappa:6.1f}: overflow/underflow")
        except Exception as e:
            print(f"  κ={kappa:6.1f}: Error - {e}")
    
    # 6.2 Edge cases
    print("\n6.2 Testing edge cases")
    test_cases = [
        (1.0001, 2.0, "a ≈ 1"),
        (1.001, 2.0, "a slightly > 1"),
        (1.01, 2.0, "a = 1.01"),
        (1.1, 2.0, "a = 1.1"),
        (10.0, 10.0, "a = 10"),
        (100.0, 100.0, "a = 100"),
        (1e6, 1e6, "a = 1e6")
    ]
    
    for a, b, desc in test_cases:
        try:
            ab = sg.otimes(a, b)
            if np.isfinite(ab) and ab > 1:
                T_ab = sg.T_transform(ab)
                T_sum = sg.T_transform(a) + sg.T_transform(b)
                error = abs(T_ab - T_sum)
                print(f"  {desc:15s}: error = {error:.2e}")
            else:
                print(f"  {desc:15s}: overflow/underflow")
        except Exception as e:
            print(f"  {desc:15s}: Error")


# =============================================================================
# 7. VISUALIZATIONS
# =============================================================================

def generate_visualizations(sg, prime_results):
    """Generate all figures for the paper"""
    print("\n" + "=" * 80)
    print("7. GENERATING VISUALIZATIONS")
    print("=" * 80)
    
    # Figure 1: Isomorphism verification
    fig1, ax1 = plt.subplots(figsize=(8, 6))
    
    a_vals = np.random.uniform(1.1, 10, 100)
    b_vals = np.random.uniform(1.1, 10, 100)
    
    left_vals = []
    right_vals = []
    
    for i in range(100):
        try:
            ab = sg.otimes(a_vals[i], b_vals[i])
            if np.isfinite(ab) and ab > 1:
                left = sg.T_transform(ab)
                right = sg.T_transform(a_vals[i]) + sg.T_transform(b_vals[i])
                if np.isfinite(left) and np.isfinite(right):
                    left_vals.append(left)
                    right_vals.append(right)
        except:
            continue
    
    ax1.scatter(right_vals, left_vals, alpha=0.6, s=30, c='blue', edgecolors='black')
    
    min_val = min(min(left_vals), min(right_vals))
    max_val = max(max(left_vals), max(right_vals))
    ax1.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Identity')
    
    ax1.set_xlabel('T(a) + T(b)', fontsize=12)
    ax1.set_ylabel('T(a ⊗ b)', fontsize=12)
    ax1.set_title(f'Isomorphism Verification (κ={sg.kappa})', fontsize=14)
    ax1.grid(True, alpha=0.3)
    ax1.legend()
    
    plt.tight_layout()
    plt.savefig('figure1_isomorphism.png', dpi=300, bbox_inches='tight')
    print("  ✓ Figure 1: Isomorphism verification saved")
    
    # Figure 2: Zeta zeros on circle
    fig2, ax2 = plt.subplots(figsize=(8, 8))
    
    radius = np.exp(1/(2*sg.kappa))
    
    theta = np.linspace(0, 2*np.pi, 100)
    ax2.plot(radius * np.cos(theta), radius * np.sin(theta), 'b-', 
             linewidth=2, label=f'|s| = e^(1/(2κ)) = {radius:.4f}')
    
    riemann_zeros_t = [14.1347, 21.0220, 25.0109, 30.4249, 32.9351]
    
    zeros_real = []
    zeros_imag = []
    for t in riemann_zeros_t:
        s = np.exp((0.5 + 1j * t) / sg.kappa)
        zeros_real.append(s.real)
        zeros_imag.append(s.imag)
    
    ax2.scatter(zeros_real, zeros_imag, color='red', s=80, marker='o', 
                edgecolors='black', label='Zeros', zorder=5)
    
    ax2.set_xlabel('Re(s)', fontsize=12)
    ax2.set_ylabel('Im(s)', fontsize=12)
    ax2.set_title(f'Zeros of ζ_κ(s) on |s| = e^(1/(2κ)) (κ={sg.kappa})', fontsize=14)
    ax2.set_aspect('equal')
    ax2.grid(True, alpha=0.3)
    ax2.legend()
    
    plt.tight_layout()
    plt.savefig('figure2_zeta_zeros.png', dpi=300, bbox_inches='tight')
    print("  ✓ Figure 2: Zeta zeros on circle saved")
    
    # Figure 3: Scale prime counting
    if prime_results and 'counting' in prime_results:
        fig3, (ax3a, ax3b) = plt.subplots(1, 2, figsize=(14, 5))
        
        x_vals = [r['x'] for r in prime_results['counting']]
        pi_vals = [r['pi_scale'] for r in prime_results['counting']]
        asym_vals = [r['asymptotic'] for r in prime_results['counting']]
        ratio_vals = [r['ratio'] for r in prime_results['counting']]
        
        ax3a.semilogx(x_vals, pi_vals, 'b-', linewidth=2, label='π_κ(x)')
        ax3a.semilogx(x_vals, asym_vals, 'r--', linewidth=2, label='T(x)/ln T(x)')
        ax3a.set_xlabel('x', fontsize=12)
        ax3a.set_ylabel('π_κ(x)', fontsize=12)
        ax3a.set_title('Scale Prime Counting Function', fontsize=14)
        ax3a.grid(True, alpha=0.3)
        ax3a.legend()
        
        ax3b.semilogx(x_vals, ratio_vals, 'g-', linewidth=2, marker='o', markersize=4)
        ax3b.axhline(y=1, color='r', linestyle='--', alpha=0.5, label='Limit = 1')
        ax3b.set_xlabel('x', fontsize=12)
        ax3b.set_ylabel('π_κ(x) / (T(x)/ln T(x))', fontsize=12)
        ax3b.set_title('Convergence to Prime Number Theorem', fontsize=14)
        ax3b.grid(True, alpha=0.3)
        ax3b.legend()
        
        plt.tight_layout()
        plt.savefig('figure3_scale_primes.png', dpi=300, bbox_inches='tight')
        print("  ✓ Figure 3: Scale prime counting saved")
    
    # Figure 4: Zeta function values
    fig4, ax4 = plt.subplots(figsize=(8, 6))
    
    s_plot = np.linspace(1.1, 3.0, 20)
    zeta_vals = []
    for s in s_plot:
        val = sg.zeta_scale(s)
        if np.isfinite(val):
            zeta_vals.append(val)
        else:
            zeta_vals.append(float('nan'))
    
    ax4.plot(s_plot, zeta_vals, 'bo-', linewidth=2, markersize=6)
    ax4.set_xlabel('s', fontsize=12)
    ax4.set_ylabel('ζ_κ(s)', fontsize=12)
    ax4.set_title(f'Scale Zeta Function ζ_κ(s) (κ={sg.kappa})', fontsize=14)
    ax4.grid(True, alpha=0.3)
    ax4.set_ylim([0.99, 1.01])
    
    plt.tight_layout()
    plt.savefig('figure4_zeta_function.png', dpi=300, bbox_inches='tight')
    print("  ✓ Figure 4: Zeta function saved")
    
    plt.show()
    print("\n  ✓ All figures saved successfully")


# =============================================================================
# MAIN EXECUTION
# =============================================================================

def main():
    """Main execution function"""
    print("\n" + "=" * 80)
    print("SCALE GROUP: COMPLETE NUMERICAL VERIFICATION SUITE")
    print("=" * 80)
    
    # Initialize scale group
    sg = ScaleGroup(kappa=config.KAPPA)
    print(f"\nκ = {sg.kappa}")
    print(f"Identity element e = {sg.identity:.6f}")
    
    # Run all verifications
    group_results = verify_group_axioms(sg)
    iso_results = verify_isomorphism(sg)
    power_results = verify_powers_and_roots(sg)
    zeta_results = verify_zeta_functions(sg)
    prime_results = verify_scale_primes(sg)
    verify_additional_tests(sg)
    
    # Generate visualizations
    generate_visualizations(sg, prime_results)
    
    # Final summary
    print("\n" + "=" * 80)
    print("FINAL VERIFICATION SUMMARY")
    print("=" * 80)
    
    print("\n✓ GROUP AXIOMS:")
    print(f"  • Associativity mean error: {group_results['associativity']['mean']:.2e}")
    print(f"  • Commutativity mean error: {group_results['commutativity']['mean']:.2e}")
    print(f"  • Identity mean error: {group_results['identity']['mean']:.2e}")
    print(f"  • Inverse mean error: {group_results['inverse']['mean']:.2e}")
    
    print("\n✓ ISOMORPHISM:")
    print(f"  • Homomorphism mean error: {iso_results['homomorphism']['mean']:.2e}")
    print(f"  • Injectivity mean error: {iso_results['injectivity']['mean']:.2e}")
    print(f"  • Surjectivity mean error: {iso_results['surjectivity']['mean']:.2e}")
    print(f"  • Inverse property mean error: {iso_results['inverse']['mean']:.2e}")
    
    print("\n✓ POWERS AND ROOTS:")
    print(f"  • Powers (n=2) mean error: {power_results['powers']['n=2']['mean']:.2e}")
    print(f"  • Roots (n=2) mean error: {power_results['roots']['n=2']['mean']:.2e}")
    
    print("\n✓ ZETA FUNCTIONS:")
    print(f"  • ζ_κ(s) = ζ(κ ln s) mean error: {zeta_results['zeta_identity']['mean']:.2e}")
    print(f"  • Zero circle mean error: {zeta_results['zeta_zeros']['mean']:.2e}")
    
    print("\n✓ SCALE PRIMES:")
    print(f"  • Prime correspondence mean error: {prime_results['prime_errors']['mean']:.2e}")
    
    # Overall assessment
    all_passed = (
        group_results['associativity']['mean'] < config.PRECISION and
        group_results['commutativity']['mean'] < config.PRECISION and
        iso_results['homomorphism']['mean'] < config.PRECISION and
        power_results['powers']['n=2']['mean'] < config.PRECISION and
        power_results['roots']['n=2']['mean'] < config.PRECISION and
        zeta_results['zeta_zeros']['mean'] < config.PRECISION
    )
    
    print("\n" + "=" * 80)
    if all_passed:
        print("✓ ALL TESTS PASSED SUCCESSFULLY - CODE READY FOR SUBMISSION")
    else:
        print("⚠ SOME TESTS FAILED - PLEASE CHECK THE RESULTS")
    print("=" * 80)
    
    return {
        'group': group_results,
        'isomorphism': iso_results,
        'powers': power_results,
        'zeta': zeta_results,
        'primes': prime_results
    }


if __name__ == "__main__":
    results = main()